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Abstract 

Tight binding (TB) models are one approach to the quantum mechanical many particle 
problem. An important role in TB models is played by hopping and overlap matrix elements 
between the orbitals on two atoms, which of course depend on the relative positions of the 
atoms involved. This dependence can be expressed with the help of Slater-Koster parameters, 
which are usually taken from tables. Recently, a way to generate these tables automatically 
was published. If TB approaches are applied to simulations of the dynamics of a system, also 
derivatives of matrix elements can appear. In this work we give general expressions for first 
and second derivatives of such matrix elements. Implemented in a computer program they 
obviate the need to type all the required derivatives of all occuring matrix elements by hand. 

1 Introduction 

Tight binding is an approach to the quantum mechanical many particle problem particularly 
valuable in cases where exact analytic solutions are not available (i.e. almost always) and other 
numerical approaches like CI or DFT are too time-consuming. For an overview see e.g. [1,2]. The 
general structure of the TB equations can be derived from DFT [1,3,4]. As is not uncommon 
in quantum mechanics, matrix elements of operators and overlap matrix elements between states 
occur also in TB. The distinctive feature of TB is that these matrix elements are considered to 
arise from aioToic- like orbitals. The precise way in which these matrix elements are obtained 
defines the type of TB approach. The matrix elements can be actually calculated from orbitals 
localised at atoms, with the orbitals obtained in some way before. Another possibility is to consider 
the matrix elements as disposable parameters (empirical TB). The matrix elements depend on the 
relative positions of the atoms involved. Regardless whether the matrix elements are obtained from 
atom-centred orbitals or introduced as parameters, it is important that this dependence on the 
relative position has the basic characteristics brought about by the properties of atomic orbitals. 
In particular, the angular dependence of the overlap of two orbitals located at different atoms is 
determined by the angular momentum quantum numbers of the orbitals. For two-centre matrix 
elements, i.e. those which depend on the positions of two atoms only, this angular dependence 
can be expressed in terms of Slater-Koster coefficients [5], which are published in various tables, 
e.g. [5,6]. Many TB schemes actually restrict to two-centre integrals, neglecting matrix elements 
with three or more centres. For a discussion of this see for example [1]. With increasing angular 
momentum quantum number I the length of the corresponding Slater-Koster tables and also the 
length of individual entries in such tables increase rapidly. Therefore a procedure to calculate the 
Slater-Koster coefficients automatically is very useful. Such an approach has been presented in [7]. 
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Figure 1: The geometric situation. Atom 1 is located at the origin. The quantities K^,Ry,R^ are 
the Cartesian coordinates of atom 2. 

If the TB approach is apphed to studies of the dynamics of a system (e.g. a molecule), also 
derivatives of the matrix elements with respect to the positions of the atoms occur, more precisely 
derivatives of first and second order [8] . To an even higher degree than for the matrix elements 
themselves, the expressions to be handled (typed into a computer) quickly become numerous and 
very complicated with increasing I. Therefore general expressions for these derivatives are most 
helpful. In this work, building on [7], we present such expressions, which, while they may still look 
a bit awkward, can be well implemented in a TB code. 

In section|2we summarise results of [7] and establish the starting point for the subsequent sections 

01 and ^ which contain our results for the first and second derivatives, respectively, of a general 
two-centre atomic-like matrix element. Some problems with the coordinates chosen are discussed 
in section |31 

2 Matrix elements 

For the convenience of the reader and also to introduce our notation this section summarises results 
of [7] on the general form of Slater-Koster matrix elements. 

Let us consider two atoms, 1 at position ri, 2 at position r2, and a quantum state on each of the 
atoms, characterised by angular and magnetic quantum numbers h,mi and ^2,7712, respectively. 
These quantum numbers determine the symmetry properties of the states relevant to the subse- 
quent discussion, whereas other quantum numbers play no explicit role. Let ri2 :— f2 — n. —■ 
1^12!^ —: be the connecting vector of atom 1 and 2. In the two-centre approximation, all matrix 
elements appearing in a tight binding approach depend only on the positions of two atoms and can 
be rephrased to depend only on the connecting vector. Therefore one of the atoms can be taken 
to be situated in the origin of the coordinate system, as illustrated in figure ^ The Cartesian 
coordinates of atom 2 then are and a, (3 are the Euler angles of the rotation bringing 

the z-axis into alignment with the connecting vector, a is defined as the rotation angle about 
the z-axis starting from the positive x-axis, < a < 27r; /3 gives the rotation angle about the 
new y'-axis (obtained as result of the a- rotation), starting from the positive z-axis, < /3 < tt. 
Note that our convention for the first Euler angle differs from [7]. A possible rotation about the 
connecting vector by a third Euler angle is irrelevant for our purposes and will not be considered. 
Between the Euler angles and the Cartesian coordinates the following relations hold (with R = 



2 



COS/3 sin/3 =-^1-^-^ (1) 



7?^ K>y 
cosa = ^= sina = ^= if {R^f + {Ryf^Q (2) 

The case (i?^)^ + (i?^)^ implies /3 = or /? = tt; a is undefined in this case and we wiU look 
at this problem in sectional 

The central point is to express a spherical harmonic Yim{Q, ^) of the unrotated coordinates (i.e. 6 is 
measured from the positive z-axis, ip from the positive x-axis) as a linear combination of spherical 
harmonics f') in the rotated frame {9' from positive ^-axis, ip' from the new x-axis ): 

fie', := Y,^ieie', ^'), ^{e', ^')) = E ^Ti^iB (3) 

Only functions of the same I are involved, because the operator effecting the rotation by the Euler 
angles 

i?(a,/3,7) = e^'^'='^e^^"'^e^'^^" (4) 

(where y denotes the intermediate y-direction after the a-rotation, z' — the final z-direction) 
commutes with [9]. From Q we have 



{yiBm^) = e™'^ {YlB\e^-'^^YiJ^ =: e™'^d^„, (/3)e^™" (5) 
where dl^^,{p) is the Wigner d-function. For our purposes 7 can always be chosen equal to 0, so 

m' 

Like [7] we use the phase convention Yi*^^{9, ip) = (— l)™y;_m(6', ip) and define real- valued spherical 
harmonics 



Y 
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Y, 



IQ 



V2(-l)™Rern,„| (7) 



l\m\ 

= \/2(-l)"Tmll|™| 



m > : Yi,, 
m < : Yir. 

Relations ^ can be summarized as 

Yirn - SmoYio + (1 - (S™o)%/2(-l)'" [r(TO)Rey,|„.| + r(-TO)Imy,|,^|] (8) 

T(m) is the discrete Heaviside function 

^<'")={; i;;::^: (9) 

In the rotated coordinate system we consider "fundamental" matrix elements of a spherically 
symmetric operator H (this includes H=I for simple overlaps between wavefunctions) 

YlfBi0^,p^)f^{\f\)Rf2i\f- R^\)Yl^J^^{02,^2)d\ = ihl2\mi\)Sm,m, (10) 

where (/i?2|'7ii|) is defined by pU|l . These quantities satisfy 

(;i/2|™i|) = (-l)'^"'^(/2Zi|mi|) (11) 
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Here the exchange ^ I2 refers to an interchange of atoms and imphes ^ — s- — ^. Because of the 
relation Hll|l not aU the matrix elements are independent. It is therefore sufficient to specify the 
elements for which li < I2', the number of independent parameters in a tight binding computer 
program is thus reduced. 

The Slater-Koster tables express a general matrix element between unrotated functions 

J Yi,rmfiB-f2Yi,m,d\ ^ {hmi\K\l2m2) (12) 

in terms of the fundamental matrix elements H1U|) and coefhcients depending on the angular mo- 
mentum quantum numbers of the states involved. The general form of the relation between ()10|l 
and H12|l according to [7] is 



{hmi\E\l2m2){a,P,R) = ^ 

m' — '. 

■ ihhlm'DiR) + 2A™,(a)A™,(a)cifj„^|o(/3)(if^^|o(/3)Gi/20)(i?) 



(13) 



where 



f (— 1)™ rT(m) cos(|m|Qf) + t(— rn) sindmla)] if to 7^ 
A„(a):=| L iim = (1^) 

„ , , f (— 1)™ [t(— to) cos(|TO|a) — r(TO) sindTOla)] if m 7^ _^ 

BM--^[ ^ ifTO = o 

^Lm' •= ^m[(— 1)" +'^{m|-m'] (1^) 

^mm' ■= Brn{l ^ ^mo) [(^l)™ C^{m|m' ~ '^Iml-m'] i^'^) 

An explicit form of the Wigner d-function is 

d'_,(/3) = 2-'(-l)'-'"' [{I + m)\{l - m)\{l + m')\{l - to')!] ^ 

(-1)'=(1 -cos/3)'-'=-'^(l + cos/3)'=+'H^ (18) 



E 

k—ky. 



k\(l — m — ky.{l — to' — fc)!(TO + m' + fc)! 



with = min(Zi, Z2), ^< = min(Z — 171,1 ~ m'), /c> ~ max(0, —to — m'). Expression p8|) has been 
obtained from [10], replacing (cos(/3/2))2 with (1 + cos/3)/2 and (sin(/3/2))2 with (1 - cos/3)/2. 
The summation limits in (|18(l are obtained from imposing the existence condition on the factorials. 
They differ from [7] but they assure the minimum number of terms in the sum. As is shown in [11] 
these conditions guarantee that the function is always defined. 

We recover the classical notation of fundamental matrix elements replacing the quantum numbers 
I and TO with their spectroscopic notation {I — {0, 1, 2, . . .} {s,p, d, . . .}, \m\ = {0, 1,2,...} 
{cr, TT, (5, . . .}). For the evaluation of cos(mck) or sin(mQ:) in Am and B„i the following relations 
(Moivre) are useful: 



cos(TOa) = E(~l)'' ( 9I) (^^" 



k=0 



-2k 

(19) 



sin(TOa)= ^ i-irf j" ){smar+'{cosay 

k=0 \ ' + / 



-2k-l 



Here [x] denotes the largest integer not larger than x. 
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3 First derivative of matrix elements 

The matrix elements have been expressed in (|13|l as functions of a, (3, R. In tight binding molec- 
ular dynamics simulations the derivatives of the matrix elements with respect to the Cartesian 
components of the connecting vector are required. We express a, f3, R as functions of R^, R^, R^ . 
As P = arccos(i?^/i?) we have 



df3 



1 



S^" - {R^R°-)I{R^) 



OR" - {R-lRf y R i?3 J ^R2 _ (_R^)2 

and, if i?^ ^ 0, tana ~ R^ /R^ yields a ~ arctan(i?^/i?^) + ip, where 



(20) 



R"" >o,Ry >0 

if^ { 27r R'^ >0,Ry <0 



da 

'dR" 



'r- 



(21) 



i + {Ry/R^Y \R^ [R'^f \ {R^Y + {RyY r'^-{r^Y 

The case i?^ = 0, i?^ 7^ implies cos a = 0, meaning a — n /2 ot a — 'dir /2] the last two expressions 
in PHI are valid in this case, too. 

If i?^ = and Ry = the connecting vector is aligned along the z-axis, and R^ — ±R. Conse- 
quently a is undefined and there are also potential problems in H2U|1 . The situation at R^ — ±i? is 
considered in section [S] The expressions in the present section are valid everywhere except at the 
poles R' = ±R. 

Let us introduce the notation F{a, /3, R) := {limi \ H \l2'm2) {a, (3, R). The derivative then is 



d 
dR^ 



(ZiTOil H \l2m2) 



OF da dF 8(3 OF OR 
da OR" ^ 'dp^m^ ^ 'dRdR^ ~ 



(i?^)2 -h {Ryf da ^i?2 - (i?^)2 9/3 R dR 



-171^ (22) 



From the derivatives of F are 



da~ ^ 



m' — 1 



ds: 



mim' QI2 



da 



ds: 



dT 



da 



da 



pt 



dT 



m2m' 



da 



A A 

dt._Jhl20) + 2A„,^^^^d'^ 



+ 2^A„.rf;j„,|od;^,|o(^i^20) + 2A„, 



E 



ds: 



h 

m2m' 



jJq, |mi|0 Imal 

dT'^ , 

mim' pl2 



Ti 



da "1"' da da "i"' da 

2|mi|i?„,,A„,dj;^^l„dj,^„^lo(/iZ20) + 2|TO2|A„ii?™2<,,|o<2|o('i^20) 



{hl2\m'\) 



(23) 



dF 
'dp 



E 

m' — 1 



'^m2m' ' '~'mim,' Qp 



m2m' 



dp 



m2Tn' ' 7nim' 



nil 



m2m' 



dp 



{hl2\m'\) 



2A„^^Am2 ^p 



'"^'"dfl ,,{hl20) + 2A„,,A^2d'^^^^,^f^{hl20) 



(24) 



dF 
dR 



E 



0/1 QI2 I T^^l 

mim' m2m' - — ' " 



mim' m2m' 



d(^i/2|m'|) 
dR 



OA A rl' 



d(;i?20) 



|mi|0 ImalO ^p 



(25) 
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The derivatives of S'^^/ and T^^^, with respect to a,P are given by 



mm 

da 



da 



(-1)™ d^\jn\m' + d\m\~m' 



= \m\B„ 



(-l)™ d^\rn\m' + '^jm|- 



55i 



dp 

dTl^^i dB„ 



An 



i-iy 



i' \'m\7n \m\—m 



d/3 



d/3 



1 — <5mO 



9q! da 

= - \m\Am ( 1 - ^rnO 



(-1)™ dimim/ - 



(-1)™ - 



(26) 
(27) 

(28) 



where use has been made of 



Bm 1 — <5mO 



(-1)' 



,dd 



\m\m' 



dd 



\m\- 



d(3 



d/3 



dAr, 



\m\Br, 



dBn 



da da 
The derivatives of the Wigner c?-function can be expressed as 



m\A^ 



ddl 



d/3 



{I + m'){l - m! + I) 



1 



{l-m'){l + m' + 1) 



"■mm' + l 



(29) 



(30) 



(31) 



and are defined for \m!\ < I. The radial part of our function is given by the fundamental matrix 
elements defined in p(J|l . which are model dependent quantities. Their derivatives with respect to 
R are thus also model dependent. 



4 Second derivative of matrix elements 

Relations (I30() and H31|l are recursive relations for computing the derivatives of Am, Bm, c^mm'i 
so in principle we can compute any higher derivative of the matrix elements. The tight binding 
approaches we are aware of require derivatives with respect to the nuclear positions up to second 
order. Below we list the corresponding expressions for these second order derivatives, valid at all 
points except the poles, which will be discussed in section [S] 



dR^dR" ^^'^^^ ' ^""^^ ~ dR^dR" da ^ ydR^ da^ dR^ dl3da ~^ R dRda)^ 
d^/3 dF dp f da d^F dp d^F R'' d^F \ 
dRf'dR" 'dp ^ dR^ [d^ dadp ^ d^'dp^ ^ 'R dRdp J ^ 



R R3 J dR ' R \dR'' dadR ' dR'' dpdR ^ R dR^J ^^^^ 



where the derivatives involved are given by 

d^P _ S^'^R'^R^ + S'^^R^R^ ~ 2R'-R''R^ {R'R'' - R^S'") {R'' - R'6' 

dR^~ RS/W^Wf ^i?2_(^.)2]i^2 

d'^a S'^'^Sy - 5y*'5''"' {R^'Sy - RyS""") + RvSv'') 



(33) 
(34) 
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aR2 



OA A rl' 



mi^m2"|TOi|0 



E 

m'=l 



{hhW\) (35) 



do? 



2\mim2\BmiBm2 - {ml + ma)^^! A 



m' = l 



m2 

i2 



|mi|0'^|m2|o(^1^20) 

ihhlm'l) (36) 



—94 A 

+ E 

m'=l 



oil \c'2 , oil 

^^2 mim' j m2m' ' mim' \Q^2^m2m' 



S^l.^, 1+21 —S]^ 



\d(3 



m\m' 



QP 7772771' 



^^2 7771777' j 7772'W' » 



nil f .^i-h ^ io/^-^t'i \f—T 

mim'\Qp2 7772777' I I J \ 8 f3 ' 



i2 

777 2 777' 



(Zi/2|m'|) 



(37) 



=2 



TOl 1-8^,^7772 



I m2 1^7771-87, 



d^«|777i|0^'*|7772|0 + ^t] 777i |0 '^l 7772 1 



+ 



+ 



'< r 

E 

m' = l 



\ dad 13 "^"^ 



I _|_ c'l [ c'2 

' ' 7772 777' ^^^'l 7772 777' 



— 

7(11777' 



(/1/2O) 
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9 



"i"'Ma/3 "2"7 V^a^/? """7 



( (9^ 



(/iZ2|m'|) 



(38) 



dRda 



m2|A„i-B^2 ]c«j^77i|0^|m2|0d^('l^20)+ ^ 



^ oil \ <-.(2 

; 777 2 777' 



+ 5"i ( g'a 

^ 777 1 777' I 7772771' 



8 \ / 8 

Q 7717777' j 7772 777' 7771777' I Q 7712777' 



di? 



(;ii2|m'| 



(39) 



8^F 



^1 m 1 1 0) '^jm2 1 + '^jmi 1 ^jm2 1 



^(W20) 



+ E 
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_|_ T^U I " T~'^2 



dR 



(40) 



7 



(— 1)™ d|„j|,^, + dj^i 



(— 1)™ + d{„j|_ 



(-1) 



m' \m\m' 



2^/ 



I m I — m ' 



d/32 



9q!2 da2 



1 " SmQ 



(-1)™ d}^jn\m' " d-lml-r 



■m^B,n [ 1 - 5„M 



(-1)™ - rf{,„|_y 



a2ri 



9/32 



,dV 



i9ai9/3 da 



\m\Br, 



(-ir 



|m|m' 



|m|- 



d/32 



d/32 



,ddl 



|m|r 



drf 



|m| — r 



d/3 



d/3 



(-ir 



,ddl 



|m|m' 

d/3 



ddl 



i 

|m| — m' 

"d^S 



(41) 
(42) 

(43) 
(44) 

(45) 



mm' 

dfida 



m,,v 

dadp 



dBm 
da 



(-1) 



,dd[ 



m' \m\m' 



dd' 



i 

\m\—m' 



\m\Ar, 



->mO 



i-iy 



d/3 
,dd\ 



d/3 



\m\7n' 



(46) 



d/3 



d/3 



5 The poles R = ±R' 

If the second atom is located on the z-axis, the azimuthal angle a cannot be defined. The derivative 
of a with respect to or diverges as the connecting vector approaches the z-axis, i.e if 
\R^\ —> R. In the same limit dp/dR^ — + and the derivatives of (3 with respect to R^ or R^ 
remain bounded, but do not converge. Related behaviour is found for the second derivatives as 
well. These problems are of course due to a well known coordinate singularity, but, given that we 
need some expressions at the poles (to type into a computer, eventually), we have to look at this 
case. A further issue to be addressed is how to handle the dependence of F and of its derivatives 
on a at the poles, i.e. what value of a to plug into the expressions where this angle is not defined. 
In order to solve both problems it is helpful to recall what exactly the Euler angles mean, see 
section 12 If both atoms are on the z-axis, a rotation about this axis does not change the matrix 
element between any two orbitals, because both orbitals get rotated by the same angle. Therefore 
there is no a-dependence for the matrix element, as can also be verified by evaluating H13I) in the 
limit /3 — > 0,7r. A partial derivative is the limit of a differential quotient. At the poles this means 
the difference between the value of the matrix element at some point P and the value at the pole, 
divided by the difference in the respective Cartesian coordinate we are considering, taken in the 
limit that P approaches the pole. But, because our point of reference is a pole, the choice of P, 
even if at an "infinitesimal" distance from the pole, determines the angle a by which to rotate 
around the z-axis, before a rotation around the thus obtained j/'-axis by the angle /3. Therefore, 
regardless along which direction in the xy-plane we move away from the pole, we are always in the 
situation shown in figure |2J this is also evident from the azimuthal symmetry about the pole. For 
/3 we thus have in this case /3 = arccos(i?^/i?) with R = \J {R^y + {R^Y, and it follows 

dR^ ~ dRy ~ dR^ ^ ^i?2 _ (i?^)2 i?2 ^B^^ r ^ > 
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Figure 2: Geometry for the derivatives at the poles. For d/dR^ the x- and z- axis are as shown, 
the y-axis, as is also indicated, is pointing towards the reader. For d/dK" a rotation about the 
z-axis by 7r/2 is necessary, and instead of the x- and y-axis we obtain the x'- and j/'-axis, the latter 
pointing towards the reader. As far as /3 is concerned, both situations are identical. The reader 
should note that is always nonnegative by construction. 




where the limit is for /3 — > 0, tt, respectively. The result 

(48) 

is as expected. Note that the total derivative of (3 exists everywhere, but is not continuous at the 
poles; thus there will be no total second derivative of (3 at the poles. 

The meaning of the Euler angles leads to the following choice for the values of a to be used in 
partial first derivatives of F at the poles: F and dF/dR do not depend on a at the poles, as follows 
from (|13|l . so its choice is arbitrary. For expressions involving dF/d/3, the choice to be made is 
a ~ for a;-derivatives and a — tt/2 for y-derivatives. The same choices of a, applied to H2U|) for 
8(3/ dR^ and dp/dR^, respectively, reproduce the result (|47ll . 

Any deviation from the pole can be entirely expressed in terms of changes in R and /3, a only enters 
to fix the direction of the deviation. Partial derivatives with respect to a, i.e. changing a while 
keeping P and R fixed, can't be made sense of at the poles. They do not occur. As furthermore 
we have dR/dR°' = ±6^°' at R^ = ±R, the general expression for the first derivatives of a matrix 
element at /3 = 0, tt respectively, is 

3 f)F 1 3F 1 f)F 

— il,,n.,lH ll.,m.) . ±6-- ± S^--^^, ± (49) 

Turning to the second partial derivatives, we find that all but one can be expressed in terms of 
derivatives of F, j3 and R. The problematic one is the mixed i?^i?^-derivative, because it involves 
two different directions in the xy-plane and thus two conflicting values of a to be used in F and 
its derivatives. All other second partial derivatives involve at most one of the variables R^ or i?^, 
and the reasoning for the first partial derivatives, as illustrated in figure HI can be carried over. 
We find at the poles 

d{R^y ' d{Rvf ' d{R^f 

d'^13 d'^13 d^(3 d^p 1 



dR^dR"^ dRvdR^' dR'^dR^ dR^dRv R^ 
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As in the case of the first derivatives, derivatives with respect to a cannot be made sense of at the 
poles. We have 

{h,mi\ H \l2,m2) = ~—j-^\a=o + -5 aaap U^O V^Oj 



' ' R'^ dl3'"-^ RdpdR 

„ , , 1 dF, 1 d^F 



{li,mi\H\l2,m2 



where the mixed derivatives are symmetric. 

For the mixed second derivative with respect to R^ and R^ we have to think of the Cartesian form 
of the expression There, from each coefficient multiplying a fundamental matrix element 

(lih]^]) we extract, if possible, a factor sinacosa(sin/3)^; this corresponds to {R^Ry)/{R?). It is 
easily seen that at the poles only those terms contribute to the derivative in question where this 
extraction is possible and where the rest of the coefficient after the extraction does not vanish. If 
we do this analysis, we find that the only nonvanishing derivatives are at /3 = 0: 

32 



(Zi,0|i7|;2,-2) - — =— + 2)(/2 + iMh - 1) {hh^)- 



dR^'dRy ' ' ' ' ' ' 2V2i?' 



(1 - 5l,^)-^^h{h + l){l2+2){l2-l) {hl2l) + 



(1 - <5,,o)(i - '^'ii)^;7f]^ + mi + i)h{h - 1) ai/22) (51) 



and 

92 



-{- 

i?2 1 2 



;V^l(^l + l)^2(/2 + l) (/l?20) - ^Mh + 1) + ;2(/2 + 1)] (?1^2l)| (52) 

At /3 = TT we obtain (— x these values. 

We think that the expressions we have given in the previous sections and completed in this one 
are suitable for implementation in a TB-based molecular dynamics code. Though the expressions 
are somewhat lengthy, they are well structured and can be coded by successive function calls. The 
functions A,„, _Bm, and df"^^, are expressed explicitly in the code, Sl^^,, T^^, and their derivatives 
with respect to a are implemented in terms of these quantities. Derivatives with respect to /3 are 
evaluated calling derivatives of the Wigner d-function, which themselves are stated as combinations 
of Wigner d-functions. 
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